$1$. A random variable $X$ is said to follow gamma$(\alpha,\beta)$ distribution.If it's probability density function given by $$\begin{equation} f(x) = \begin{cases} \frac{e^{-\frac{x}{\beta}}x^{\alpha-1}}{\beta^{\alpha}\Gamma(\alpha)}, 0 < x< \infty \\ 0~~~~~~ ,otherwise \end{cases} \end{equation}$$.
(a) . Draw $f(x)$ for different values of $\alpha$ and $\beta$. Create two plots. In one plot fix $\alpha = 2$ and vary $\beta = 2,4,6.5$. In the other plot, fix $\beta = 3$ and vary $\alpha = 4,6,7.5$. Each plot will contain three functions. Use different functions.Add appropiate legends to the plots.
(b). It is known that if $X_{1},X_{2},\cdots,X_{m} \sim$ gamma$(1,\beta)$, then $\sum_{i = 1}^m X_{i} \sim \mathcal{gamma}(m,\beta)$. Verify the statement using simulation. Use $m = 10$ and $\beta = 3$. Consider $1000$ replication to visualize the distribution of $\sum_{i =1}^m X_{i}$.
using Distributions,Plots,SpecialFunctions,QuadGK;
using LaTeXStrings;
function f(y,a,b)
(exp(-y/b)*(y^(a-1)))/((b^a)*gamma(a))
end
f (generic function with 3 methods)
a = 2
b_vals= [2,4,6.5];
y = 0:0.01:15;
p = plot()
for b in b_vals
l = [f.(y,a,b)]
p = plot!(y,l,label = "Gamma($(a),$(b))",fg_legend = false)
end
display(p)
b = 3
a_vals= [4,6,7.5]
y = 0:0.01:15;
p1 = plot()
for a in a_vals
l = [f.(y,a,b)]
p1= plot!(y,l,label = "Gamma($(a),$(b))",legend = :topleft,fg_legend = false, grid = false)
end
display(p1)
x1 = []
m = 10
beta = 3
rep = 1000
for i in 1:rep
push!(x1,rand(Gamma(m,beta)))
end
# beta density function with (m,beta) as parameters
function f(x)
(exp(-x/beta)*x^(m-1))/(gamma(m)*(beta^m))
end
f (generic function with 2 methods)
histogram(x,normalize = true,label = "histogram of Xi",c = 4,fg_legend = false)
plot!(f,10,60,label = "gamma(10,3)",c= 6, linewidth = 4,style = :dash)
$2$. Consider the following function $\begin{equation} g(x) = \begin{cases} \frac{1}{\sqrt{2x}}e^{-\frac{x^2}{4}} , 0 < x< \infty \\ 0 ~~~~~~~~, otherwise \end{cases} \end{equation}$
(a) Using Monte Carlo method compute the integral of the function on the interval $[a,b] = [2,5]$ based on $n = 1000$ samples. Call it $I_{n}$.
(b) Using inbuilt function to compute the integral, call it $I$.
(c) Print the absolute difference between $I$ and $I_{n}$.
(d) Try simulate, show that $I_{n}$ converges to $I$ as $n \to \infty$.
We want to compute the integral $\int_{2}^{5} g(x) dx$.Take $a = 2,b = 5$.Then the integral can be written as $(b-a)\int_{a}^{b}\frac{1}{b-a}\frac{1}{\sqrt{2x}}e^{-\frac{x^2}{4}}dx = E_{f}[g(X)] ,f\sim \mathcal{U}(a,b)$.
So we generate random samples $x_{i}$ from $\mathcal{Uniform(a,b))}$, then calculate $g(x_{i})$ then take the mean and multiply it by $(b-a)$. We do this process for $n = 1,2,\cdots, 1000$.
function g(x)
(1/sqrt(2*x))*(exp(-x^2/4))*(x>0)
end
g (generic function with 1 method)
a = 2; b = 5;
I_n = []
n = 1000;
for i in 1:n
x = rand(Uniform(2,5),i)
push!(I_n,(b-a)*mean(g.(x)))
end
I_n
1000-element Vector{Any}:
0.007956155850635398
0.012879270597985045
0.16372531948611774
0.10490097549397108
0.018840875047610554
0.18014317027367047
0.12268589027643673
0.2560520985104784
0.11287763607330115
0.10772582282532644
0.1291242512115906
0.1421804375174541
0.10591406591090198
⋮
0.1159685271251723
0.12555332242330125
0.12746130845837733
0.13371915904913068
0.12781254015601573
0.12663111881770747
0.1242725332912191
0.12648038568286088
0.12775120893860106
0.12413690244898444
0.12225380072351498
0.12431813357814087
mean(I_n)
0.12245044474453692
exact_integral = quadgk(g,a,b)[1]
0.12290701094708961
plot(I_n,label = "approx integral")
hline!([exact_integral],label = "exact integral")
abs_diff = []
for j in 1:length(I_n)
push!(abs_diff,I_n .- exact_integral)
end
length(abs_diff)
1000
plot([abs_diff],legend = false)
hline!([0],label =:none)